home *** CD-ROM | disk | FTP | other *** search
/ Programmer Power Tools / Programmer Power Tools.iso / microcrn / issue_40.arc / DAIMS.ARC / OCEAN.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1988-02-10  |  3.6 KB  |  156 lines
  1. #include <math.h>
  2. #include "viscosity.h"
  3. #include "matrix.hxx"
  4. #include "vimatrix.hxx"
  5. #include "Cheb_vector.hxx"
  6. #include "ekman.hxx"
  7. #include "ocpanel.hxx"
  8. #include "ocean.hxx"
  9. /* 
  10. -*++ ocean_layer::ocean_layer(): 
  11. ** 
  12. ** (*++ history: 
  13. **     16 Jan 88    Bruce Eckel    Creation date
  14. ** ++*)
  15. ** 
  16. ** (*++ detailed: 
  17. ** ++*)
  18. */
  19.  
  20. ocean_layer::ocean_layer(int initial_modes, double initial_lambda) :
  21. Qnew(initial_modes), Qold(initial_modes), h(initial_modes), 
  22.     A(initial_modes, initial_lambda)
  23. {
  24.     ocean_lambda = initial_lambda;
  25.     dt = ocean_lambda/500;
  26.     nstep = 0; 
  27.     time = 0;
  28.     ocean_nmodes = initial_modes;
  29.     A = A.inverse();
  30.     ocean_xldomain = -1.0;
  31.     ocean_xrdomain = 1.0;
  32.     ocean_viscosity_value = 0;
  33.     ocean_viscosity_type = NONE;
  34.     running = 0;
  35. }
  36.     
  37.  
  38. /* 
  39. -*++ ocean_layer::step(): step the model forward by one time increment
  40. ** 
  41. ** (*++ history: 
  42. **     16 Jan 88    Bruce &    Creation date
  43. ** ++*)
  44. ** 
  45. ** (*++ detailed: 
  46. ** ++*)
  47. */
  48. #define NL "\n"
  49. void ocean_layer::step(ekman_layer & ekman)
  50. {
  51.     if (!running) return;
  52. #if 0
  53.     if (ocean_nmodes != ekman.resolution()) {
  54.     ocean_nmodes = ekman.resolution();
  55.     nstep = 0;
  56.     time = 0;
  57.     A  = vi_matrix(ocean_nmodes,ocean_lambda);
  58.     A = A.inverse();
  59.     Qold =  Cheb_vector (ocean_nmodes);
  60.     Qnew =  Cheb_vector (ocean_nmodes);
  61.     h =  Cheb_vector (ocean_nmodes);
  62.     }
  63. #endif   
  64. #if 0
  65.     cout << "\nnstep = " << nstep << NL;
  66.     cout.flush();
  67.     cout << "h = " << h << NL;
  68.     cout.flush();
  69.     cout << "h.prime = " << h.prime() << NL;
  70.     cout.flush();
  71.     cout << "ekman.pumping_vector() = " << ekman.pumping_vector() << NL;
  72.     cout.flush();
  73.     cout << "Qold = " << Qold << NL;
  74.     cout.flush();
  75. #endif
  76.  
  77. /*    if (!nstep++) {
  78.     Qnew = Qold +  (ekman.pumping_vector() - h.prime()) * dt;
  79.     }
  80.     else 
  81.     Qnew = Qold + ( (ekman.pumping_vector() - h.prime())) * (2 * dt);
  82. */
  83.  
  84.     Qnew = Qold +  (ekman.pumping_vector() - h.prime()) * dt;
  85.     Qnew[Qnew.size() -2] = 0;
  86.     Qnew[Qnew.size() -1] = 0;
  87. //    cout << "Qnew = " << Qnew << NL;
  88. //    cout.flush();
  89.     h = A * Qnew;
  90. //    cout << "h physical" << h.physical() << NL;
  91. //    cout.flush();
  92.     Qold = Qnew;
  93. }
  94.  
  95. /* 
  96. -*++ ocean_layer::update(): get new values from the ocean panel
  97. ** 
  98. ** (*++ history: 
  99. **     16 Jan 88    Bruce &    Creation date
  100. ** ++*)
  101. ** 
  102. ** (*++ detailed: 
  103. ** ++*)
  104. */
  105.  
  106.  void ocean_layer::update(ocpanel & ocean_panel, ekman_layer & ekman)
  107. {
  108.     if (ocean_panel.nmodes() != ekman.resolution()  || ocean_panel.changed()) {
  109.     ocean_lambda = ocean_panel.lambda();
  110.     ocean_nmodes = ekman.resolution();
  111.     dt = ocean_lambda/100;
  112.     A = vi_matrix(ocean_nmodes, ocean_lambda);
  113.     A = A.inverse();
  114.     Qnew = Cheb_vector(ocean_nmodes);
  115.     Qold = Cheb_vector(ocean_nmodes);
  116.     h = Cheb_vector(ocean_nmodes);
  117.     ocean_panel.reset(); /* so it doesn't show dirty flag */
  118.     }
  119.     running = ocean_panel.run();
  120.     if (ocean_panel.restart()) {
  121.     ocean_panel.clear_restart();
  122.     ocean_nmodes = ekman.resolution();
  123.     ocean_lambda = ocean_panel.lambda();
  124.     A = vi_matrix(ocean_nmodes, ocean_lambda);
  125.     A = A.inverse();
  126.     dt = ocean_lambda/100;
  127.     nstep = 0; 
  128.     time = 0;
  129.     ocean_xldomain = -1.0;
  130.     ocean_xrdomain = 1.0;
  131.     ocean_viscosity_value = 0;
  132.     Qnew = Cheb_vector(ocean_nmodes);
  133.     Qold = Cheb_vector(ocean_nmodes);
  134.     h = Cheb_vector(ocean_nmodes);
  135.     }
  136. }
  137.  
  138.  
  139. /* 
  140. -*++ ocean_layer::diagnostics(): send results to ocpanel display
  141. ** 
  142. ** (*++ history: 
  143. **     16 Jan 88    Bruce &    Creation date
  144. ** ++*)
  145. ** 
  146. ** (*++ detailed: 
  147. ** ++*)
  148. */
  149.  
  150. void ocean_layer::diagnostics(ocpanel & ocean_panel)
  151. {
  152.     ocean_panel.display(h.physical());
  153. }
  154.  
  155.  
  156.